#include "AtlasStyle.C"
#include "AtlasUtils.h"
#include "NAssymetricError.C"

void set_jet_binlabel(TH1D *hist) {

    TAxis *x = hist->GetXaxis();
 
    for( int i=1; i <= hist->GetNbinsX(); i++) {
        ostringstream stm;
        stm << i-1;
        x->SetBinLabel(i, TString(stm.str()) );
    }
}

TCanvas* plot_mtop(TString hist_fpath, const int ijet=0, const int ierr=0) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "mtop_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    gROOT->LoadMacro("AtlasUtils.C");

    TFile *f = new TFile(hist_fpath);

    /*
    TH1D * h_data    = (TH1D*) f->Get("data10_mtop_tight"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_mtop_tight"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_mtop_tight"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_mtop_tight"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_mtop_tight"   + ht_ext); h_zjets->SetDirectory(0);
    */

    TH1D * h_data    = (TH1D*) f->Get("data10_mtop"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_mtop"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_mtop"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_mtop"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_mtop"   + ht_ext); h_zjets->SetDirectory(0);

    //TH1D * h_diboson = (TH1D*) f->Get("mc_diboson_mwt_tight_"+ ht_ext); h_diboson->SetDirectory(0);
    TH1D * h_qcd = h_data->Clone(); h_qcd->Reset();

    /*
    if (ijet == 0) {
        for (int j=1; j<4; j++) {
            TString ht_ext_tmp = "_j";
            ht_ext_tmp += j;
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_mtop_qcd" + ht_ext_tmp); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
        }
    } else {
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_mtop_qcd" + ht_ext); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
    }
    */
    TString prefix_hqcd_name = "data10_mtop_qcd";

    switch(ierr) {
        case 1:
            prefix_hqcd_name += "_estat";
            break;
      
        case 2:
            prefix_hqcd_name += "_esyst";
            break;

        default:
            break;
    }

    TH1D * h_qcd_tmp = (TH1D*) f->Get(prefix_hqcd_name + ht_ext); h_qcd_tmp->SetDirectory(0);
    h_qcd->Add(h_qcd_tmp);

    h_qcd->SetDirectory(0);

    TH1D * h_qcd_err = h_qcd->Clone(); h_qcd_err->SetDirectory(0);

    h_data->SetTitle("");
    h_data->GetXaxis()->SetTitle("MeV");
    h_data->GetYaxis()->SetTitle("Events");
    h_data->GetXaxis()->SetLabelSize(.03);
    h_data->GetYaxis()->SetLabelSize(.03);
    h_data->GetXaxis()->SetTitleSize(.03);
    h_data->GetYaxis()->SetTitleSize(.03);
    h_data->GetXaxis()->SetTitleOffset(1.5);
    h_data->GetYaxis()->SetTitleOffset(1.5);
    h_data->SetMarkerStyle(20);

    h_ttbar->SetFillColor(0);
    h_ttbar->GetXaxis()->SetTitle("MeV");
    h_ttbar->GetYaxis()->SetTitle("Events");
    h_ttbar->GetXaxis()->SetLabelSize(.03);
    h_ttbar->GetYaxis()->SetLabelSize(.03);
    h_ttbar->GetXaxis()->SetTitleSize(.03);
    h_ttbar->GetYaxis()->SetTitleSize(.03);
    h_ttbar->GetXaxis()->SetTitleOffset(1.5);
    h_ttbar->GetYaxis()->SetTitleOffset(1.5);

    h_stop->SetFillColor(kBlue-3);
    h_wjets->SetFillColor(kOrange);
    h_zjets->SetFillColor(kGreen-4);  
    //h_diboson->SetFillColor(kBlue-3);

    h_qcd->SetFillColor(kOrange+2);

    h_qcd_err->SetFillColor(kOrange+4);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();

    // adding up histogram for making stack plot
    h_ttbar->Add(h_stop);
    h_ttbar->Add(h_wjets);
    h_ttbar->Add(h_zjets);
    h_ttbar->Add(h_qcd);

    // deciding the first drawing histogram (data or MC?)
    Double_t max_data = 0;
    Double_t max_mc   = 0;

    for ( int i=1; i<=h_data->GetNbinsX(); i++) {
        Double_t cbin = h_data->GetBinContent(i);
        if ( cbin > max_data ) max_data = cbin; 
    }

    for ( int i=1; i<=h_ttbar->GetNbinsX(); i++) {
        Double_t cbin = h_ttbar->GetBinContent(i);
        if ( cbin > max_mc ) max_mc = cbin;
    }

    if ( max_data > max_mc ) {
        h_data->Draw("E1");
        h_ttbar->Draw("HIST SAME");
    } else {
        h_ttbar->Draw("HIST");
    }

    h_stop->Add(h_wjets);
    h_stop->Add(h_zjets);
    h_stop->Add(h_qcd);
    h_stop->Draw("HIST SAME");

    h_wjets->Add(h_zjets);
    h_wjets->Add(h_qcd);
    h_wjets->Draw("HIST SAME");

    h_zjets->Add(h_qcd);
    h_zjets->Draw("HIST SAME");

    h_qcd->Draw("HIST SAME");

    h_qcd_err->SetMarkerSize(0);
    h_qcd_err->SetFillStyle(3001);
    h_qcd_err->SetFillColor(kBlack);
    h_qcd_err->Draw("E2 SAME");

    h_data->Draw("E1 SAME");

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_data,"data (36.6 pb^{-1})","p");
    legend->AddEntry(h_ttbar,"t#bar{t}","f");
    legend->AddEntry(h_stop,"Single Top","f");
    legend->AddEntry(h_wjets,"W+jets","f");
    legend->AddEntry(h_zjets,"Z+jets","f");
    legend->AddEntry(h_qcd,"QCD","f");
    legend->AddEntry(h_qcd_err,"QCD uncertainty","f");
    legend->Draw();

    // KS test
    Double_t KSTest = h_data->KolmogorovTest(h_ttbar);
    cout << "KS test: " << KSTest << endl;

    TString ks_str = "KS: ";
    ks_str += KSTest;

    myText2(0.7, 0.5, 1, .03, ks_str);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    TString upflavor, uptag, upjet;
    upflavor = "[#mu+jets]";
    myText(.2,0.88, kRed+4, upflavor);

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    //f->Close();

    return c;
}

TCanvas* plot_eff_l2t_data(TString hist_fpath, const int ijet=0) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "eff_l2t_data_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    gROOT->LoadMacro("AtlasUtils.C");

    TFile *f = new TFile(hist_fpath);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();

    TH1D * h_tight = (TH1D*) f->Get("data10_d0sig_tight" + ht_ext); h_tight->SetDirectory(0);
    TH1D * h_loose = (TH1D*) f->Get("data10_d0sig_loose" + ht_ext); h_loose->SetDirectory(0);
 
    TGraphAsymmErrors *g = EfficiencyPlot(h_loose, h_tight);
    g->SetMarkerColor(kBlack);
    g->SetMarkerStyle(20);
    g->Draw("P");

    g->GetXaxis()->SetTitle("d0 significance");
    g->GetYaxis()->SetTitle("efficiency");
    g->GetXaxis()->SetLabelSize(.03);
    g->GetYaxis()->SetLabelSize(.03);
    g->GetXaxis()->SetTitleSize(.03);
    g->GetYaxis()->SetTitleSize(.03);
    g->GetXaxis()->SetTitleOffset(1.5);
    g->GetYaxis()->SetTitleOffset(1.5);

    g->SetLineColor(kBlue);
    g->SetMarkerColor(kBlue);
    g->SetMarkerStyle(22);
    g->Draw("ALP");

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(g,"data (36.6 pb^{-1})","p");

    legend->Draw();

    //myText2(0.6, 0.8, 1, .03, plot_title);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    delete f;

    return c;
}

TCanvas* plot_eff_l2t_mc(TString hist_fpath, const int ijet=0, const bool inc_mcqcd=false, const bool inc_data=false) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "eff_l2t_mc_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    gROOT->LoadMacro("AtlasUtils.C");

    TFile *f = new TFile(hist_fpath);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();

    TH1D *h_tight_mix = (TH1D*) f->Get("mc_mix_d0sig_tight" + ht_ext); h_tight_mix->SetDirectory(0);
    TH1D *h_loose_mix = (TH1D*) f->Get("mc_mix_d0sig_loose" + ht_ext); h_loose_mix->SetDirectory(0);
 
    TGraphAsymmErrors *g_mix = EfficiencyPlot(h_loose_mix, h_tight_mix);

    g_mix->GetXaxis()->SetTitle("d0 significance");
    g_mix->GetYaxis()->SetTitle("efficiency");
    g_mix->GetXaxis()->SetLabelSize(.03);
    g_mix->GetYaxis()->SetLabelSize(.03);
    g_mix->GetXaxis()->SetTitleSize(.03);
    g_mix->GetYaxis()->SetTitleSize(.03);
    g_mix->GetXaxis()->SetTitleOffset(1.5);
    g_mix->GetYaxis()->SetTitleOffset(1.5);

    g_mix->SetLineColor(kBlue);
    g_mix->SetMarkerColor(kBlue);
    g_mix->SetMarkerStyle(22);
    g_mix->Draw("ALP");

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);

    legend->AddEntry(g_mix,"MC mix","lp");

    if ( inc_mcqcd ) {

        TH1D *h_tight_qcd = (TH1D*) f->Get("mc_qcd_d0sig_tight" + ht_ext); h_tight_qcd->SetDirectory(0);
        TH1D *h_loose_qcd = (TH1D*) f->Get("mc_qcd_d0sig_loose" + ht_ext); h_loose_qcd->SetDirectory(0);
 
        TGraphAsymmErrors *g_qcd = EfficiencyPlot(h_loose_qcd, h_tight_qcd);

        g_qcd->SetLineColor(kRed);
        g_qcd->SetMarkerColor(kRed);
        g_qcd->SetMarkerStyle(21);
        g_qcd->Draw("LP");

        legend->AddEntry(g_qcd,"MC QCD (Alpgen_J1_filtmu)","lp");
    }

    if ( inc_data ) {
        TH1D * h_tight = (TH1D*) f->Get("data10_d0sig_tight" + ht_ext); h_tight->SetDirectory(0);
        TH1D * h_loose = (TH1D*) f->Get("data10_d0sig_loose" + ht_ext); h_loose->SetDirectory(0);
 
        TGraphAsymmErrors *g = EfficiencyPlot(h_loose, h_tight);
        g->SetMarkerColor(kBlack);
        g->SetMarkerStyle(20);
        g->Draw("P");

        legend->AddEntry(g,"data (36.6 pb^{-1})","p");
    }
    legend->Draw();

    //myText2(0.6, 0.8, 1, .03, plot_title);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    delete f;

    return c;
}

TCanvas* plot_nbjets(TString hist_fpath, const int ijet=0, const bool logy=false, const bool qcd=true, const int ierr=0) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "nbjets25_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    gROOT->LoadMacro("AtlasUtils.C");

    TFile *f = new TFile(hist_fpath);

    /*
    TH1D * h_data    = (TH1D*) f->Get("data10_nbjets25_tight"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_nbjets25_tight"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_nbjets25_tight"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_nbjets25_tight"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_nbjets25_tight"   + ht_ext); h_zjets->SetDirectory(0);
    */

    TH1D * h_data    = (TH1D*) f->Get("data10_nbjets25"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_nbjets25"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_nbjets25"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_nbjets25"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_nbjets25"   + ht_ext); h_zjets->SetDirectory(0);

    //TH1D * h_diboson = (TH1D*) f->Get("mc_diboson_mwt_tight_"+ ht_ext); h_diboson->SetDirectory(0);
    TH1D * h_qcd = h_data->Clone(); h_qcd->Reset();

    /*
    if (ijet == 0) {
        for (int j=1; j<4; j++) {
            TString ht_ext_tmp = "_j";
            ht_ext_tmp += j;
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_nbjets25_qcd" + ht_ext_tmp); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
        }
    } else {
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_nbjets25_qcd" + ht_ext); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
    }
    */
    TString prefix_hqcd_name = "data10_nbjets25_qcd";

    switch(ierr) {
        case 1:
            prefix_hqcd_name += "_estat";
            break;
      
        case 2:
            prefix_hqcd_name += "_esyst";
            break;

        default:
            break;
    }


    TH1D * h_qcd_tmp = (TH1D*) f->Get(prefix_hqcd_name + ht_ext); h_qcd_tmp->SetDirectory(0);
    h_qcd->Add(h_qcd_tmp);

    h_qcd->SetDirectory(0);

    TH1D * h_qcd_err = h_qcd->Clone(); h_qcd_err->SetDirectory(0);

    h_data->SetTitle("");
    h_data->GetXaxis()->SetTitle("NBjets");
    h_data->GetYaxis()->SetTitle("Events");
    h_data->GetXaxis()->SetLabelSize(.03);
    h_data->GetYaxis()->SetLabelSize(.03);
    h_data->GetXaxis()->SetTitleSize(.03);
    h_data->GetYaxis()->SetTitleSize(.03);
    h_data->GetXaxis()->SetTitleOffset(1.5);
    h_data->GetYaxis()->SetTitleOffset(1.5);
    h_data->SetMarkerStyle(20);

    h_ttbar->SetFillColor(0);
    h_ttbar->GetXaxis()->SetTitle("NBjets");
    h_ttbar->GetYaxis()->SetTitle("Events");
    h_ttbar->GetXaxis()->SetLabelSize(.03);
    h_ttbar->GetYaxis()->SetLabelSize(.03);
    h_ttbar->GetXaxis()->SetTitleSize(.03);
    h_ttbar->GetYaxis()->SetTitleSize(.03);
    h_ttbar->GetXaxis()->SetTitleOffset(1.5);
    h_ttbar->GetYaxis()->SetTitleOffset(1.5);

    h_stop->SetFillColor(kBlue-3);
    h_wjets->SetFillColor(kOrange);
    h_zjets->SetFillColor(kGreen-4);  
    //h_diboson->SetFillColor(kBlue-3);

    h_qcd->SetFillColor(kOrange+2);

    h_qcd_err->SetFillColor(kOrange+4);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();
    if (logy) c->SetLogy();

    // adding up histogram for making stack plot
    h_data->Draw("E1");

    h_ttbar->Add(h_stop);
    h_ttbar->Add(h_wjets);
    h_ttbar->Add(h_zjets);
    if (qcd) h_ttbar->Add(h_qcd);
    h_ttbar->Draw("HIST SAME");

    h_stop->Add(h_wjets);
    h_stop->Add(h_zjets);
    if (qcd) h_stop->Add(h_qcd);
    h_stop->Draw("HIST SAME");

    h_wjets->Add(h_zjets);
    if (qcd) h_wjets->Add(h_qcd);
    h_wjets->Draw("HIST SAME");

    if (qcd) h_zjets->Add(h_qcd);
    h_zjets->Draw("HIST SAME");

    if (qcd) h_qcd->Draw("HIST SAME");

    if (qcd) {
        h_qcd_err->SetMarkerSize(0);
        h_qcd_err->SetFillStyle(3001);
        h_qcd_err->SetFillColor(kBlack);
        h_qcd_err->Draw("E2 SAME");
    }

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_data,"data (36.6 pb^{-1})","p");
    legend->AddEntry(h_ttbar,"t#bar{t}","f");
    legend->AddEntry(h_stop,"Single Top","f");
    legend->AddEntry(h_wjets,"W+jets","f");
    legend->AddEntry(h_zjets,"Z+jets","f");
    if (qcd) legend->AddEntry(h_qcd,"QCD","f");
    if (qcd) legend->AddEntry(h_qcd_err,"QCD uncertainty","f");
    legend->Draw();

    // KS test
    //Double_t KSTest = h_data->KolmogorovTest(h_ttbar);
    //cout << "KS test: " << KSTest << endl;

    //TString ks_str = "KS: ";
    //ks_str += KSTest;

    //myText2(0.7, 0.5, 1, .03, ks_str);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    TString upflavor, uptag, upjet;
    upflavor = "[#mu+jets]";
    myText(.2,0.88, kRed+4, upflavor);

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    // showing each bin value data:qcd
    for (Int_t i=0; i < h_data->GetNbinsX() + 1; i++) {
        Double_t N_sel     = h_data->GetBinContent(i);
        Double_t N_qcd     = h_qcd->GetBinContent(i);
        Double_t N_qcd_err = h_qcd->GetBinError(i);
        cout << i << ":" << N_sel << " " << N_qcd << " +- " << N_qcd_err << endl;
    }

    delete f;

    return c;
}

TCanvas* plot_njets(TString hist_fpath, const int ijet=0, const bool logy=false, const bool qcd=true, const int ierr=0) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "njets25_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    gROOT->LoadMacro("AtlasUtils.C");

    TFile *f = new TFile(hist_fpath);

    /*
    TH1D * h_data    = (TH1D*) f->Get("data10_njets_tight"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_njets_tight"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_njets_tight"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_njets_tight"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_njets_tight"   + ht_ext); h_zjets->SetDirectory(0);
    */

    TH1D * h_data    = (TH1D*) f->Get("data10_njets"       + ht_ext); h_data->SetDirectory(0); set_jet_binlabel(h_data);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_njets"   + ht_ext); h_ttbar->SetDirectory(0);set_jet_binlabel(h_ttbar);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_njets"    + ht_ext); h_stop->SetDirectory(0); set_jet_binlabel(h_stop);
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_njets"   + ht_ext); h_wjets->SetDirectory(0);set_jet_binlabel(h_wjets);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_njets"   + ht_ext); h_zjets->SetDirectory(0);set_jet_binlabel(h_zjets);

    TH1D * h_qcd = h_data->Clone(); h_qcd->Reset();set_jet_binlabel(h_qcd);

    /*
    if (ijet == 0) {
        for (int j=1; j<4; j++) {
            TString ht_ext_tmp = "_j";
            ht_ext_tmp += j;
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_njets_qcd" + ht_ext_tmp); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
        }
    } else {
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_njets_qcd" + ht_ext); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
    }
    */

    TString prefix_hqcd_name = "data10_njets_qcd";

    switch(ierr) {
        case 1:
            prefix_hqcd_name += "_estat";
            break;
      
        case 2:
            prefix_hqcd_name += "_esyst";
            break;

        default:
            break;
    }

    TH1D * h_qcd_tmp = (TH1D*) f->Get(prefix_hqcd_name + ht_ext); h_qcd_tmp->SetDirectory(0);
    h_qcd->Add(h_qcd_tmp);

    h_qcd->SetDirectory(0);

    TH1D * h_qcd_err = h_qcd->Clone(); h_qcd_err->SetDirectory(0); set_jet_binlabel(h_qcd_err);

    h_data->SetTitle("");
    h_data->GetXaxis()->SetTitle("Njets");
    h_data->GetYaxis()->SetTitle("Events");
    h_data->GetXaxis()->SetLabelSize(.03);
    h_data->GetYaxis()->SetLabelSize(.03);
    h_data->GetXaxis()->SetTitleSize(.03);
    h_data->GetYaxis()->SetTitleSize(.03);
    h_data->GetXaxis()->SetTitleOffset(1.5);
    h_data->GetYaxis()->SetTitleOffset(1.5);
    h_data->SetMarkerStyle(20);

    h_ttbar->SetFillColor(0);
    h_ttbar->GetXaxis()->SetTitle("Njets");
    h_ttbar->GetYaxis()->SetTitle("Events");
    h_ttbar->GetXaxis()->SetLabelSize(.03);
    h_ttbar->GetYaxis()->SetLabelSize(.03);
    h_ttbar->GetXaxis()->SetTitleSize(.03);
    h_ttbar->GetYaxis()->SetTitleSize(.03);
    h_ttbar->GetXaxis()->SetTitleOffset(1.5);
    h_ttbar->GetYaxis()->SetTitleOffset(1.5);

    h_stop->SetFillColor(kBlue-3);
    h_wjets->SetFillColor(kOrange);
    h_zjets->SetFillColor(kGreen-4);  
    //h_diboson->SetFillColor(kBlue-3);

    h_qcd->SetFillColor(kOrange+2);

    h_qcd_err->SetFillColor(kOrange+4);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();
    if (logy) c->SetLogy();

    // adding up histogram for making stack plot
    h_data->Draw("E1");

    h_ttbar->Add(h_stop);
    h_ttbar->Add(h_wjets);
    h_ttbar->Add(h_zjets);
    if (qcd) h_ttbar->Add(h_qcd);
    h_ttbar->Draw("HIST SAME");

    h_stop->Add(h_wjets);
    h_stop->Add(h_zjets);
    if (qcd) h_stop->Add(h_qcd);
    h_stop->Draw("HIST SAME");

    h_wjets->Add(h_zjets);
    if (qcd) h_wjets->Add(h_qcd);
    h_wjets->Draw("HIST SAME");

    if (qcd) h_zjets->Add(h_qcd);
    h_zjets->Draw("HIST SAME");

    if (qcd) h_qcd->Draw("HIST SAME");

    if (qcd) {
        h_qcd_err->SetMarkerSize(0);
        h_qcd_err->SetFillStyle(3001);
        h_qcd_err->SetFillColor(kBlack);
        h_qcd_err->Draw("E2 SAME");
    }

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_data,"data (36.6 pb^{-1})","p");
    legend->AddEntry(h_ttbar,"t#bar{t}","f");
    legend->AddEntry(h_stop,"Single Top","f");
    legend->AddEntry(h_wjets,"W+jets","f");
    legend->AddEntry(h_zjets,"Z+jets","f");
    if (qcd) legend->AddEntry(h_qcd,"QCD","f");
    if (qcd) legend->AddEntry(h_qcd_err,"QCD uncertainty","f");
    legend->Draw();

    // KS test
    //Double_t KSTest = h_data->KolmogorovTest(h_ttbar);
    //cout << "KS test: " << KSTest << endl;

    //TString ks_str = "KS: ";
    //ks_str += KSTest;

    //myText2(0.7, 0.5, 1, .03, ks_str);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    TString upflavor, uptag, upjet;
    upflavor = "[#mu+jets]";
    myText(.2,0.88, kRed+4, upflavor);

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    // showing each bin value data:qcd
    for (Int_t i=0; i < h_data->GetNbinsX() + 1; i++) {
        Double_t N_sel     = h_data->GetBinContent(i);
        Double_t N_qcd     = h_qcd->GetBinContent(i);
        Double_t N_qcd_err = h_qcd->GetBinError(i);
        cout << i << ":" << N_sel << " " << N_qcd << " +- " << N_qcd_err << endl;
    }

    delete f;

    return c;
}

TCanvas* plot_mwt(TString hist_fpath, const int ijet=0, const bool qcd=true, const bool mc_qcd=false, const double mc_qcd_scale=1.0, const int ierr=0) {

    gROOT->SetStyle("Plain");       //no grey background on plots
    gStyle->SetOptStat(0);
    SetAtlasStyle();

    TObjArray *finfo = hist_fpath.Tokenize("/");

    TString ht_ext = "_j";
    ht_ext += ijet;

    TString hist_fname = ((TObjString *) finfo->At( finfo->GetEntries()-1 ))->GetString();

    TString plot_title = "mwt_" + hist_fname + ht_ext;
    plot_title.ReplaceAll(".root","");

    gROOT->LoadMacro("AtlasUtils.C");


    ostringstream st_scale_mc2d;
    st_scale_mc2d << mc_qcd_scale;

    TFile *f = new TFile(hist_fpath);

    /*
    TH1D * h_data    = (TH1D*) f->Get("data10_mwt_tight"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_mwt_tight"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_mwt_tight"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_mwt_tight"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_mwt_tight"   + ht_ext); h_zjets->SetDirectory(0);
    */

    TH1D * h_data    = (TH1D*) f->Get("data10_mwt"       + ht_ext); h_data->SetDirectory(0);
    TH1D * h_ttbar   = (TH1D*) f->Get("mc_ttbar_mwt"   + ht_ext); h_ttbar->SetDirectory(0);
    TH1D * h_stop    = (TH1D*) f->Get("mc_stop_mwt"    + ht_ext); h_stop->SetDirectory(0); 
    TH1D * h_wjets   = (TH1D*) f->Get("mc_wjets_mwt"   + ht_ext); h_wjets->SetDirectory(0);
    TH1D * h_zjets   = (TH1D*) f->Get("mc_zjets_mwt"   + ht_ext); h_zjets->SetDirectory(0);

    TH1D * h_mcqcd   = h_data->Clone(); h_mcqcd->Reset();
    TH1D * h_mcqcdbb = h_data->Clone(); h_mcqcdbb->Reset();

    if ( mc_qcd ) {
        h_mcqcd   = (TH1D*) f->Get("mc_qcd_mwt"     + ht_ext); h_mcqcd->SetDirectory(0);
        h_mcqcdbb = (TH1D*) f->Get("mc_qcdbb_mwt"   + ht_ext); h_mcqcdbb->SetDirectory(0);

        h_mcqcd->Scale(mc_qcd_scale);
        h_mcqcdbb->Scale(mc_qcd_scale);
    }

    //TH1D * h_diboson = (TH1D*) f->Get("mc_diboson_mwt_tight_"+ ht_ext); h_diboson->SetDirectory(0);
    TH1D * h_qcd = h_data->Clone(); h_qcd->Reset();

    /*
    if (ijet == 0) {
        for (int j=1; j<4; j++) {
            TString ht_ext_tmp = "_j";
            ht_ext_tmp += j;
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_mwt_qcd" + ht_ext_tmp); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
        }
    } else {
            TH1D * h_qcd_tmp = (TH1D*) f->Get("data10_mwt_qcd" + ht_ext); h_qcd_tmp->SetDirectory(0);
            h_qcd->Add(h_qcd_tmp);
    }
    */

    TString prefix_hqcd_name = "data10_mwt_qcd";

    switch(ierr) {
        case 1:
            prefix_hqcd_name += "_estat";
            break;
      
        case 2:
            prefix_hqcd_name += "_esyst";
            break;

        default:
            break;
    }

    TH1D * h_qcd_tmp = (TH1D*) f->Get(prefix_hqcd_name + ht_ext); h_qcd_tmp->SetDirectory(0);
    h_qcd->Add(h_qcd_tmp);

    h_qcd->SetDirectory(0);

    TH1D * h_qcd_err = h_qcd->Clone(); h_qcd_err->SetDirectory(0);

    h_data->SetTitle("");
    h_data->GetXaxis()->SetTitle("MeV");
    h_data->GetYaxis()->SetTitle("Events");
    h_data->GetXaxis()->SetLabelSize(.03);
    h_data->GetYaxis()->SetLabelSize(.03);
    h_data->GetXaxis()->SetTitleSize(.03);
    h_data->GetYaxis()->SetTitleSize(.03);
    h_data->GetXaxis()->SetTitleOffset(1.5);
    h_data->GetYaxis()->SetTitleOffset(1.5);
    h_data->SetMarkerStyle(20);

    h_ttbar->SetFillColor(0);
    h_ttbar->GetXaxis()->SetTitle("MeV");
    h_ttbar->GetYaxis()->SetTitle("Events");
    h_ttbar->GetXaxis()->SetLabelSize(.03);
    h_ttbar->GetYaxis()->SetLabelSize(.03);
    h_ttbar->GetXaxis()->SetTitleSize(.03);
    h_ttbar->GetYaxis()->SetTitleSize(.03);
    h_ttbar->GetXaxis()->SetTitleOffset(1.5);
    h_ttbar->GetYaxis()->SetTitleOffset(1.5);

    h_stop->SetFillColor(kBlue-3);
    h_wjets->SetFillColor(kOrange);
    h_zjets->SetFillColor(kGreen-4);  
    //h_diboson->SetFillColor(kBlue-3);

    h_qcd->SetFillColor(kOrange+2);
    h_qcd_err->SetFillColor(kOrange+4);

    h_mcqcd->SetFillColor(kOrange+2);
    h_mcqcdbb->SetFillColor(kOrange+4);

    TCanvas * c = new TCanvas(plot_title,plot_title,1);
    c->SetTickx();
    c->SetTicky();

    // adding up histogram for making stack plot

    h_ttbar->Add(h_stop);
    h_ttbar->Add(h_wjets);
    h_ttbar->Add(h_zjets);
    if (qcd)    h_ttbar->Add(h_qcd);
    if (mc_qcd) {
        h_ttbar->Add(h_mcqcd);
        h_ttbar->Add(h_mcqcdbb);
    }
    // deciding the first drawing histogram (data or MC?)
    Double_t max_data = 0;
    Double_t max_mc   = 0;

    for ( int i=1; i<=h_data->GetNbinsX(); i++) {
        Double_t cbin = h_data->GetBinContent(i);
        if ( cbin > max_data ) max_data = cbin; 
    }

    for ( int i=1; i<=h_ttbar->GetNbinsX(); i++) {
        Double_t cbin = h_ttbar->GetBinContent(i);
        if ( cbin > max_mc ) max_mc = cbin;
    }

    if ( max_data > max_mc ) {
        h_data->Draw("E1");
        h_ttbar->Draw("HIST SAME");
    } else {
        h_ttbar->Draw("HIST");
    }

    h_stop->Add(h_wjets);
    h_stop->Add(h_zjets);
    if (qcd) h_stop->Add(h_qcd);
    if (mc_qcd) {
        h_stop->Add(h_mcqcd);
        h_stop->Add(h_mcqcdbb);
    }
    h_stop->Draw("HIST SAME");

    h_wjets->Add(h_zjets);
    if (qcd) h_wjets->Add(h_qcd);
    if (mc_qcd) {
        h_wjets->Add(h_mcqcd);
        h_wjets->Add(h_mcqcdbb);
    }
    h_wjets->Draw("HIST SAME");

    if (qcd) h_zjets->Add(h_qcd);
    if (mc_qcd) {
        h_zjets->Add(h_mcqcd);
        h_zjets->Add(h_mcqcdbb);
    }
    h_zjets->Draw("HIST SAME");

    if (qcd) h_qcd->Draw("HIST SAME");

    if (qcd) {
        h_qcd_err->SetMarkerSize(0);
        h_qcd_err->SetFillStyle(3001);
        h_qcd_err->SetFillColor(kBlack);
        h_qcd_err->Draw("E2 SAME");
    }

    if (mc_qcd) {
        h_mcqcd->Draw("HIST SAME");
        h_mcqcdbb->Draw("HIST SAME");
    }

    h_data->Draw("E1 SAME");

    legend = new TLegend(0.7,0.57,0.95,0.85);
    legend->SetName("legend");
    legend->SetHeader(plot_title);
    legend->SetFillColor(0);
    legend->SetFillStyle(0);
    legend->SetLineColor(0);
    legend->SetBorderSize(0);
    legend->SetTextSize(.03);
    legend->AddEntry(h_data,"data (36.6 pb^{-1})","p");
    legend->AddEntry(h_ttbar,"t#bar{t}","f");
    legend->AddEntry(h_stop,"Single Top","f");
    legend->AddEntry(h_wjets,"W+jets","f");
    legend->AddEntry(h_zjets,"Z+jets","f");
    if (qcd) legend->AddEntry(h_qcd,"QCD","f");
    if (qcd) legend->AddEntry(h_qcd_err,"QCD uncertainty","f");
    if (mc_qcd) {
        legend->AddEntry(h_mcqcd,"QCD*" + TString(st_scale_mc2d.str()),"f");
        legend->AddEntry(h_mcqcdbb,"QCDbb*" + TString(st_scale_mc2d.str()),"f");
    }
    legend->Draw();

    // KS test
    Double_t KSTest = h_data->KolmogorovTest(h_ttbar);
    cout << "KS test: " << KSTest << endl;

    TString ks_str = "KS: ";
    ks_str += KSTest;

    myText2(0.7, 0.5, 1, .03, ks_str);

    //ATLAS_LABEL(0.6,0.88);
    //myText( 0.72,0.88,1,"Preliminary");

    TString upflavor, uptag, upjet;
    upflavor = "[#mu+jets]";
    myText(.2,0.88, kRed+4, upflavor);

    gPad->RedrawAxis();
    //gPad->Range(0., 87., 0., 160000.);

    c->SaveAs(plot_title + ".eps");
    c->SaveAs(plot_title + ".png");
    c->SaveAs(plot_title + ".pdf");

    delete f;

    return c;
}
